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Abstract. The existence of two stationary solutions of the nonlinear Boltzmann 
equation for inelastic hard spheres or disks is investigated. They are restricted neither 
to weak dissipation nor to small gradients. The one-particle distribution functions 
are assumed to have an scaling property, namely that all the position dependence 
occurs through the density and the temperature. At the macroscopic level, the state 
corresponding to both is characterized by uniform pressure, no mass flow, and a linear 
temperature profile. Moreover, the state exhibits two peculiar features. First, there is 
a relationship between the inelasticity of collisions, the pressure, and the temperature 
gradient. Second, the heat flux can be expressed as being linear in the temperature 
gradient, i.e. a Fourier-like law is obeyed. One of the solutions is singular in the 
elastic limit. The theoretical predictions following from the other one are compared 
with molecular dynamics simulation results and a good agreement is obtained in the 
parameter region in which the Fourier state can be actually observed in the simulations, 
namely not too strong inelasticity. 
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1. Introduction 

The knowledge of dilute granular gases has improved greatly in the last years [H El E]. 
Among many others, one of the fundamental reasons for studying granular gases is that 
they are considered as a proving ground for kinetic theory and non-equilibrium statistical 
mechanics. The underlying idea is that the grains can be assimilated to atoms or 
molecules and, therefore, granular gases to molecular gases. Due to the inherent energy 
dissipation in collisions between grains, there is no equilibrium state for granular gases 
and any state of them is of a non-equilibrium nature. At a phenomenological level, 
granular gases exhibit many similarities with ordinary gases, although also some strong 
differences jl]. 

Grains are often modeled as inelastic hard spheres or disks and, in the simplest 
versions, tangential friction is neglected. For this kind of models, the inelastic Boltzmann 
kinetic equation [5J provides the accurate starting point to study low-density granular 
gases. The validity of this equation is justified on exactly the same grounds as its 
elastic, molecular, limit. It can be derived either heuristically or from the Liouville 
equation in the asymptotic small density limit [6] . Moreover, the theoretical predictions 
derived from it have been found to be in very good agreement with molecular dynamics 
simulation results, when the comparison is carried out inside the appropriate range of 
the parameters defining the system. The above includes the velocity distribution of the 
simplest possible state of a granular gas, the so-called homogeneous cooling state [TJ, 
and its exponential tail [HI E], the inelastic hydrodynamic equations to Navier-Stokes 
order, and the expressions for the transport coefficients appearing in them JTUJ ITT] . 
More severe conditions such as the initial departure of a granular gas from homogeneity 
due to an instability have also been investigated |12| . finding again that the Boltzmann 
equation accurately predicts the behavior of the system in the low density limit. 

The Chapman-Enskog procedure for deriving hydrodynamic equations, assumes the 
existence of a special kind of solutions of the Boltzmann equation, the normal solutions, 
and that they can be constructed by means of an expansion in a formal inhomogeneity 
parameter [T3l [Ti] . The hydrodynamic Navier-Stokes equations correspond to keeping 
up to the first order in the gradients contributions to the distribution function, giving rise 
to second order in the gradients terms in the hydrodynamic equations, when employed in 
the formal expressions of the heat and momentum fluxes. A peculiarity of granular gases 
is that hydrodynamic gradients are often induced by inelasticity. This holds particularly 
in the case of steady states, since the only way of compensating for the energy dissipation 
in collisions is through energy fluxes associated to gradients in the system. Consequently, 
in these systems there is a coupling between gradients and inelasticity, so that the 
restriction to small gradients, as it is the case of the Navier-Stokes equations, implies 
in many cases the restriction to small inelasticity too. For some states, like the steady 
shear flow, the situation is even more complex, since they are inherently rheological and 
the Navier-Stokes approximation never applies [To] . 

Of course, both limitations mentioned above, the assumption of a normal solution 
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and of an expansion in powers of the non-uniformity of the system, are overcome if the 
exact solution of the Boltzmann equation for a given state is known. In this case, all 
the transport properties can be computed with no assumption and/or approximation. 
The problem is that very little is known about exact inhomogeneous solutions of the 
Boltzmann equation that are relevant for transport. For molecular, elastic systems, the 
available exact solutions correspond to the very unrealistic case of Maxwell molecules 
[T6l [PT] . On the other hand, a few years ago the existence of an inhomogeneous exact 
solution of the Boltzmann equation for smooth inelastic hard spheres or disks was 
suggested [IB] . It has the property that all the spatial dependence in the velocity 
distribution occurs through its second velocity moment or, equivalently, the granular 
temperature. The macroscopic state is stationary and with gradients in only one 
direction. Moreover, the relationship between the heat flow and the temperature 
gradient can always be expressed in a linear form, no matter the value of the parameters 
of the system. For this reason, it was referred to as the Fourier state. 

Here the study of the Fourier state is undertaken in more detail. At the beginning, 
the aim of this study was twofold: to construct a solution of the equations derived in [18] 
beyond the Navier-Stokes approximation considered there, and to investigate at what 
extent the predicted state can be actually observed in molecular dynamics simulations. 
The latter is a first step towards the possibility of seeing the state in real experiments. 
Nevertheless, the analysis to be presented indicates the possible existence of two different 
solutions of the inelastic Boltzmann equation having the properties associated to the 
Fourier state as described above. The first of them agree, in the appropriate range 
of parameters, with the solution obtained in the Navier-Stokes approximation following 
from the Chapman-Enskog solution to the (inelastic) Boltzmann equation. On the other 
hand, the other solution can not be inferred from a Chapman-Enskog-like algorithm. 
The reason is that it is singular, in the sense that its elastic limit does not correspond to 
any solution of the elastic Boltzmann equation. For similar reason, this solution is not 
captured either by the hydrodynamic Navier-Stokes equations for dilute granular gases. 

The remainder of the paper is organized as follows. The Fourier state of a dilute 
granular gas is defined in Sec. where the assumed form of the one-particle distribution 
function for this state is reviewed. Substitution of its expression into the inelastic 
Boltzmann equation and scaling of the velocity, leads to a closed kinetic equation 
in which all the explicit space dependence has been eliminated. In Sec. [3] a Sonine 
expansion of the distribution function solution of the Boltzmann equation is considered. 
The expansion is truncated to the lowest order leading to nontrivial contributions to 
the even and odd parts of the distribution function. Equations for the coefficients of the 
remaining terms are obtained by taking velocity moments in the Boltzmann equation. 
In these equations only the lowest nonlinear contributions are kept. Then, two different 
solutions, referred to as "regular" and "singular" , respectively, are identified. The origin 
and consistency of the singular solution is discussed in Sec. HI The purpose there is not 
to prove rigourously the existence of the solution, but to check the consistency of the 
calculations and approximations indicating that this is the case. 
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A comparison of the theoretical predictions with molecular dynamics simulation 
results is carried out in Sec. [5j The state described by the singular solution of the 
Boltzmann equation has never been observed. A possible and probable reason for it 
is that this solution corresponds to a highly unstable state. On the other hand, for 
weak inelasticity, a good agreement is found in the bulk of the simulated systems with 
the predictions from the regular solution. As the inelasticity is increased, the Sonine 
approximation fails and later on, the Fourier state becomes very difficult to reach, in 
a significant space region of the system, in the simulations. Section [6] presents a brief 
discussion of the results, possible extensions, and also some open questions. 



2. The Fourier state 



Consider a dilute granular gas composed of smooth inelastic hard spheres (d = 3) or 
disks (d = 2) of diameter a and mass m, being a the velocity-independent coefficient of 
normal restitution. The one-particle distribution function f(r,v,t) is assumed to obey 
the inelastic Boltzmann equation [5]. For steady distributions depending on position 
only through the x-coordinate, it has the form 

v x ^f(x,v) = J[v\f], (1) 

where the Boltzmann collision term J is a functional of the form 

J[v\f] = a d - X J dvt J da9(g-a)g-a[a- 2 f(x,v')f(x,v[) 

-f(x,v)f(x,v 1 )}. (2) 

In the above expression, g = v — V\ is the relative velocity, <r is a unit vector along the 
line of centers of the two colliding particles at contact, 6 is the Heaviside step function, 
and v', v[ are the precollisional velocities leading to v and V\. Moreover, attention here 
will be restricted to solutions having a vanishing average velocity, i.e., 

dv vf(x, v) = 0. (3) 

Balance equations are obtained in the usual way by taking velocity moments in Eq. (TjQ), 

S^ = 0, (4) 

+ T(x)C(x) = 0, (5) 
where n(x) and T(x) are the number density and the temperature, respectively, 

n(x) = [dv f{x, v), ^ (x f (x) = \ dv^ f(x, v), (6) 



and the pressure tensor and heat flux q are defined as 

/[ THJlP' 
dvmViVjf(x,v), q= dv——vf(x,v). (7) 
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Finally, ((x) is the cooling rate due to the energy dissipation in collisions. Its expression 
is 



rrnr 2 o d 1 (l — a 2 ) 



dv \ dv x g f(x,v)f(x,vi). 



(8) 



AdT (^) nT 

Solving Eq. ([1]) for given boundary conditions, is a very hard task. In ref . [18] , the 
search of a solution having the scaling form 



n[x) 



m 



2T(x) 



d/2 



(p(c) 



ni 



2T(x) 



1/2 



(9) 



was proposed. This expression is similar to the distribution function of the so-called 
homogeneous cooling state [5], but with the position x playing the role of the time. 
The above functional form had already been used in ref. [IS] to fit the molecular 
dynamics results for the distribution function of the steady state of a granular gas with 
a temperature gradient. Equations (Tj0) and ([9]) imply that the pressure tensor and, 
therefore, the pressure p = ^iPu/d = nT are uniform. Due to Eqs. (J4]) and ([6]), the 
function (p must verify 

J dcp(c) = 1, J dccp(c) =0, J dcc 2 p(c) = -. (10) 

Use of Eq. ([9]) into Eq. ([5]), taking into account Eqs. (0) and (E]), yields 
dT 



with 



dx 



pa 



d-i 



2\ d ~ 1 

a )7r 2 



J dc j dci \c — ci\ 3 (p(c)ip(ci) 



(11) 
(12) 

(13) 

Moreover, Eq. ffTTj) establishes a relationship between the pressure, the temperature 
gradient, and the coefficient of restitution a (through the function ip), 

9 ' ^ (14) 



2r(^) Jdcc 2 c x ip(c) 
It follows that the temperature profile is strictly linear in x, 
dT n 

—— = a = constant, 
dx 



d-l 



pa 

The heat flux in the x direction is 





~2T{x)~ 


1/2 , 


q x = 




"J 




m 





/ 2T X 1/2 J dcc 2 Cx ^ dT 



p dec c x ip(c) = I — 



-T- (15) 



\m J a d ~ l I[(p\ dx 
This has the form of a Fourier law, with the heat flux coupled linearly to the temperature 
gradient. It is worth to stress that here it has been derived without any explicit 
restriction to small gradients. This is the reason why this state was referred to as 
the Fourier state [15] . 

To get a closed equation for the function ip, Eq. is substituted into the 

Boltzmann equation (Q, taking into account Eq. (TTTT) . The result reads: 



c 

I[p] { C x p{c) + y— • [C<p(c)] 



a 



,1-d 



j[cM. 



(16) 
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From a mathematical point of view, the problem is fully analogous to the identification 
of the distribution function of the homogeneous cooling state [5]. 



3. Sonine expansion 

Because of symmetry, the function (p(c) must be an even function of the vector 
component of c perpendicular to the x-axis, c±. Then, (p is expanded in series of 
Sonine polynomials as 



oo oo 



i=0 j=0 



(17) 



Here, the function has been decomposed into its even and odd in c x parts, and the Sn^ 
are the Sonine polynomials [13], [H] , 

:-i) fc i>+* + i) , 



k=0 



T(n + k + l)(i- k)\k\ 



They verify the orthogonality relation 

dzz n e-*S!f)(z)SU\z) 



r(n + z + l) 



(19) 



io *i 
where T(z) is the Gamma function and dij the Kronecker delta. In particular, it is 

S®(z) = l, S<»(z)=n + l-z, (20) 

for all n. Using the relation ffl9|) . it is obtained that the necessary and sufficient 
conditions for ip to verify Eqs. ( flOl) are 

a 00 = 1, b 00 = 0, aio + (d- l)a 01 = 0. (21) 

In the following, the expansion given in Eq. ( fTTl) will be considered to the lowest 
nontrivial order, keeping only the terms proportional to a 01 (and a w ), b i, and b 10 , 
and neglecting the remaining terms. This corresponds to the so-called first Sonine 
approximation [TJ]. Explicitly, the approximation considered can be expressed as 

'd 



ip(c) 



1 - a i(c - dc x ) + 



1 1 3 L 

— boi + - Oio ) <■,■ 



b ic x c 2 - (&io - boi)C 



(22) 



It is worth to emphasize that this approximation does not imply by itself the isotropy 
of the state. For instance, it is 



2 1 (d - l)a i 
dc c 2 xV (c) = - + 



dc c 2 ± (p(c) 



1 



(1 - a i) , (23) 



indicating that the diagonal elements of the pressure tensor are not the same if aoi does 
not vanish. Also, the dimensionless heat flux 



Qx 



m 


1/2 , 

Si- 


_2T(x)_ 


p J 



(24) 
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can be decomposed into the two components 

Q xx = [ dcc 3 Mc) = -^, (25) 



Q x ± = J dcc x c 2 Mc) = - {d 4 1)&Q1 . (26) 

This implies the anisotropy in the energy being carried in the x direction. To identify 
the coefficients aoi, &oi> and b\o, velocity moments are taken in Eq. (fl6l) . using the 
expansion of <p(c) in the first Sonine approximation. The equation contains terms that 
are nonlinear in tp(c), leading consequently to contributions that are nonlinear in the 
coefficients to be determined. Here the simplification is made of keeping only terms 
up to second degree in the coefficients a i, &oi> an d &io, i- e - neglecting those terms 
proportional to aoi&oi^io wr ^h n i + ra 2 + "-3 > 3. A similar approximation to determine 
the coefficient characterizing the first Sonine approximation to the distribution function 
of the homogeneous cooling state [5j [7] , was found to lead to a very good estimate, even 
for strong inelasticity [21 [20]. In the present case, the accuracy of the approximation 
being used is just assumed for the sake of simplicity, without further justification. 

Consider I[(p] defined in Eq. f fl2l . A quadratic calculation as indicated above gives 

J dcj dc 1 \c-c 1 \McMc 1 )= y > [1 + 3(601,610)], (27) 

n(h h , _ 456f + 18(d - 1)601610 + 3(d 2 - 1)6^ 

9{h ^ ho) = 64d(d + 2)(d + 4) (28) 

and then, using Eqs. f[24|) - f[26l) . it follows that 



(1 - « 2 W^2 5 / 2 

= rW2) [(d-i) 6oi + 3M 11 + 610)1 ' (29) 

Now, Eq. fflBT) is multiplied by c 2 and integrated over c. After some algebra, it is 
obtained 

, x (d- l)(2d + 3-3a) 
(1 - a)(d - l)(6 i - 36i ) = - a 01 ± >- 

x [(d - l)6oi + 36i ] • (30) 

Upon obtaining this expression, Eq. ( |28i) has been employed. Notice that an independent 
relationship between the coefficients can not be derived by multiplying Eq. ([TEI) by c\ 
and posterior integration over c. The reason is that the balance equation for the energy 
has already been employed in the derivation of Eq. (1161) . Therefore, in order to get 
two more independent equations for the coefficients aoi, 601 , and 6i , higher velocity 
moments of Eq. ( iTTl) have been considered. In the Appendix, the equations obtained 
by multiplying the kinetic equation by c\ and by c 2 c x , and integrating afterwards over 
the velocity c, are reported. Then the parameters a i, 601, and 6i can be obtained by 
numerically solving the system formed by Eqs. f[30l) . flA.lj) . and flA.2j) . The equations 



are of second degree and have two different sets of real solutions. They will be referred 
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Figure 1. The dimensionless parameters aoi, &io, and &oi determining the distribution 
function of the "regular" Fourier state in the first Sonine approximation considered 
in the text, as a function of the coefficient of normal restitution a. The solid lines 
correspond to d — 2 and the dashed ones to d = 3. 

to as "regular" solution and "singular" solution, respectively. The reason for this 
nomenclature will be discussed later on. 

In Fig. [1] the parameters a i, &oi> an d b w are plotted as a function of a for the 
regular solution. The results corresponding to both d = 2 and d = 3 are shown. The 
parameters defining the singular solution are given in Fig. [2j It is seen that aoi does 
not vanish in any of the two solutions for a ^ 1, then indicating the anisotropy of the 
pressure tensor, as discussed above. In the limit a — > 1, there is a strong difference 
between the behaviors of the regular and the singular solutions. While in the former the 
three coefficients tend to vanish, in the latter they tend to finite, non-vanishing values. 
Moreover, for the regular solution a m becomes very small and 6 i — ^io> i n the limit of 
small inelasticity (a close to one). This is consistent with the results obtained to Navier- 
Stokes order by the Chapman- Enskog procedure in the first Sonine approximation, that 
leads to a distribution function verifying a Q i = 0, b i = b w [21] . When the inelasticity 
of the system increases, the values of &oi and bio for the regular solution grow both 
very fast and the difference between them becomes significant. On the other hand, the 
behavior of the singular solution for a — > 1 strongly differs from the Chapman-Enskog 
solution, since not only the three coefficients tend to non-zero values, but the limits of 
&oi and &io are definitely different, having even opposite sign. A direct consequence of 
this behavior is that, while the first solution tends to the elastic equilibrium Gaussian 
for a — > 1, the other one does not. This is the reason to qualify them as regular and 
singular solutions, respectively. 
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Figure 2. The dimensionless parameters aoi, 610, and &oi determining the distribution 
function of the "singular" Fourier state in the first Sonine approximation considered 
in the text, as a function of the coefficient of normal restitution a. The solid lines 
correspond to d — 2 and the dashed ones to d = 3. 



4. The singular solution 

It has been mentioned above that a peculiarity of the singular solution of the Boltzmann 
equation for the Fourier state is that the coefficients 6qi an d b w do not vanish in 
the elastic limit, tending to values with opposite signs in the elastic limit. A direct 
implication of this, following from Eqs. (125]) and ff26l) is that the two components of the 
heat flux, Q xx and Q x ±, also have opposite signs. That means that particles moving 
in one direction have on the average larger values of v x , while particles moving in the 
opposite direction have, also on the average, larger values of v\. Although one can be 
prompted to conclude that something fundamental is being violated by this solution, 
we have not been able to identify any argument leading to discard it a priori. A trivial 
first test of the consistency of the singular solution is that the net heat flux must vanish 
in the elastic limit, since the dissipation and, therefore, the temperature gradient vanish 
in it. In Fig. [3] the dimensionless heat flux Q x corresponding to the singular solution is 
plotted as a function of the coefficient of normal restitution. It is seen to vanish when 
a —>■ 1, for both d = 2 and d — 3, as it should. 

In order to get additional information about whether the singular solution is an 
artifact of the Sonine approximation carried out in sec. [3j it seems worth to investigate 
the asymptotic behavior of the solutions of the Boltzmann equation (Tl6l) in the limit 
a — > 1. Assume that in this limit, the asymptotic behavior of the solutions of the 
equation reads 

^( C )~^°)( C ) + eV {1 )(c), (31) 
where e = (1 — a 2 ) 1 / 2 , and the parameter q > and the functions ip^°'(c) and ^ l \c) 
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Figure 3. The dimensionless heat flux Q x for the singular solution of the inelastic 
Boltzmann equation for the Fourier state. The solid line correspond to d = 2 and the 
dashed one to d = 3. 



are to be determined in the following by consistency. As a consequence, 

Q a ~ QW + <*QW , (32) 



with 



Similarly, 



fl^> = /*^<c). (33) 



d-1 

7M = 2I T^ /dc/ddlc-diVcMcO-T^ + eV^, (34) 



d-l 

7 (°) = ^ E 3y /dc/ddlc-dlV^^d), (35) 



7 (i) 



7T 2 



y rf c y ic - ciiv o) ( c )p (i) ( c i)- (36) 



r(f) 

The asymptotic behavior of the inelastic Boltzmann collision operator is easily identified 

as 

J[c\g}~ jM[c\g} + e 2 jW[c\g], (37) 

where is the elastic Boltzmann collision operator and the explicit for of [c\g] 
will not be needed in the following. 

Putting e = in Eq. (fTBj) it is obtained 

Qf J<°W>]=0. (38) 

Therefore, either 

Qf = (39) 
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or 

J(°)[c|^°)] =0. (40) 

Since Eq. ( 1401) implies that <^°)(c) is the Gaussian for which the dimensionless heat 
flux vanishes, it is concluded that Eq. (139]) is always verified. When also Eq. (140]) 
is fulfilled, the behavior corresponding to the regular solution, and to the Chapman- 
Enskog expansion, is recovered. The other possibility, namely [c\(p^] ^ but still 
Qx°^ = 0, is expected to lead to the singular solution. To investigate further this issue, 
the next order balance in the Boltzmann equation ( |T6l) is considered. It is given by 

e W°) j^°)(c) + \^ ■ [epV>(c)] } * eV 1 -^ J(°)[c|^(°)]. (41) 

If Qx^ 7^ 0, balance of this equation requires that q = 2. Then multiplication of the 
equation by c 2 x and integration over c yields 

^-^-i^Vl (42) 

with 

«2 = /^ (43) 

Similarly, multiplication by Cj_ and integration over c gives 

7 (o )Q (o) = _ 2g a) ff i--y dcc^c^ *], (44) 

where = — QaS- Summation of Eqs. (1411 and ( 1441 . taking into account that 
the collision operator conserves the kinetic energy, leads to Eq. fl39|) providing a 
consistency test for the existence of the singular solution. Moreover, since the singular 
solution differs from the gaussian, there is no reason to expect that the right hand side 
of Eqs. ( 142]) and (l44t) vanish. In summary, the asymptotic behavior of the solutions of 
the Boltzmann equation is qualitatively consistent with the existence of the "singular" 
solution discussed above. Going to higher order terms trying to actually determine the 
asymptotic values of the moments of the singular solutions seems a formidable task. 



5. Molecular Dynamics simulation results 

The Fourier state described in the previous sections is a "bulk state" , in the sense that 
it is expected to be shown by a dilute granular gas far away from the boundaries under 
the appropriate conditions, e.g. stationarity, no macroscopic mass flow, and gradients in 
only one direction [191 122] . Whether or not this scenario can be generated with enough 
accuracy in experiments with real boundary conditions, or even in particle simulations 
with idealized ones, is not at all trivial. In principle, it could be thought that the bulk 
state would be reached by increasing the size of the system, while suitably scaling its 
properties [23]. Nevertheless, due to the coupling between the gradients, the pressure, 
and the inelasticity in the Fourier state, the identification of the right way of scaling 
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the system turns out to be far from trivial. Although several different attempts were 
made trying to reach the Fourier state corresponding to the singular solution in the 
simulations, we did not succeed, possible because it is highly unstable. Consequently, 
attention will be restricted in the following to the regular solution, and all the references 
to theoretical predictions in the remaining of this section must be understood as dealing 
with that solution, although no established explicitly for the sake of brevity. 

The macroscopic one-dimensional state we are considering is known to exhibit a 
hydrodynamic instability, leading to the development of transversal inhomogeneities 
[2U [25]. In a two-dimensional system with periodic boundary conditions in the y 
direction, the instability is controlled by the aspect ratio of the system L y /L x , where 
L x and L y are the dimensions of the system in the direction of the gradients and 
perpendicular to it, respectively. For given values of all the other parameters, the 
instability shows up when the aspect ratio exceeds a given critical value. Of course, it 
is possible that other instabilities occur in different regions of the space defined by the 
physical parameters of the system. 

To investigate the validity of the theoretical predictions presented in the previous 
sections and its accessibility, molecular dynamics (MD) simulations of a system of 
inelastic hard disks (d = 2) have been performed. The simulations started with 
the particles uniformly distributed on a square lattice and with a Gaussian velocity 
distribution. Several ways of injecting energy through the walls were investigated, and 
it was concluded that the most efficient one to generate a bulk region in the system, 
are the so-called thermal walls [261 ETj . In these walls, particles are absorbed at the 
surface and instantaneously reemitted with a velocity determined by an equilibrium 
Maxwell-Boltzmann distribution, with a given wall temperature. The simulation data 
reported in the following have been obtained with two thermal walls located at x = 
and x = L x , respectively, and with periodic boundary conditions in the y direction. 

In all the cases to be reported, it was found that, after a transient period, the 
system reached a macroscopic steady state with gradients only in the x direction and no 
flow field. The results presented below have been time averaged, once the system was 
in the steady state, and also over several independent trajectories. The temperature 
parameters of the thermal walls were fixed as follows. First, the temperature To of the 
wall located at x = was arbitrarily fixed. Then, the temperature of the wall at x = L x 
was determined by using the theoretical prediction for the temperature profile, Eq. (flT|) . 
In this calculation, it was taken into account that the temperature of the thermal wall 
at x = (x = L x ) is expected to correspond to the temperature of the gas extrapolated 
to x = — A(0) (x = L x + X(L X )), where X(x) is the local mean free path. 

Figure 0] shows the pressure, temperature, and density profiles measured in a system 
with a = 0.99. In agreement with the theoretical predictions, it is observed that the 
pressure is uniform and the temperature profile lineal, outside a boundary layer. In this 
case, the bulk of the system, arbitrarily identified as the region in which the temperature 
profile is linear, extends over most of it. A much more demanding and fundamental check 
of the theory is presented in Fig. [51 There, the marginal velocity distribution <p x (c x ) 



The Fourier state of dilute granular gases 



13 



defined by 

ip x (c x ) = J dc y (p{c) (45) 

is plotted for three different values of x, namely for the layers centered at x — 275<r, 
475a, and 675<r, respectively. The results for the other layers in the bulk are similar. 
Of course, in each layer the velocities have been scaled with the local temperature. 
It is seen that the overlap of the data is very good over a wide range of velocities, 
hence confirming the scaling of the distribution function in Eq. ([9]). Moreover, the 
difference with a Gaussian (dotted line) is significant, while the theoretical prediction 
derived (solid line) here provides a much more accurate description of the distribution 
in the thermal velocity region, roughly \c x \ < 2 . To investigate the tails of the velocity 
distribution, the latter is also plotted on a logarithmic scale. As expected, the theoretical 
prediction derived here fails to correctly describe the tails, since it is based on the first 
Sonine approximation. On the other hand, the scaling seems to be verified within the 
statistical uncertainties. 

Finally, in Fig. [6l the dimensionless heat flux defined in Eq. (|24"|) is plotted as a 
function of the position of the layer considered. Although for this property the statistical 
error, defined as the standard deviation of the measured values, is rather large, a quite 
uniform value along the system is neatly observed. Moreover, this value is in good 
agreement with the theoretical prediction following from Eqs. (1251) and (1261) . with the 
coefficients determined by Eqs. fl30l . flA.lj) . and flA.2j) . 

Similar results have been obtained for other values of the coefficient of restitution in 
the interval 0.9 < a < 1, although the bulk region in which the Fourier state shows up 
becomes narrower as a decreases. In addition, the first Sonine approximation becomes 
less accurate. Nevertheless, there is strong evidence that the assumed Fourier state 
is present. To illustrate these comments, in Fig. [7J the temperature and heat flux 
profiles are plotted for a system with a = 0.9. Now the bulk of the system, identified 
as the region in which the temperature profile is linear and the heat flux uniform, is 
restricted to an interval roughly between 250cx and 550cr. The reduced distribution 
function tp(c) for three different positions inside the bulk region is shown in Fig. [HJ 
Significant discrepancies with the theoretical prediction in the first Sonine approximation 
are observed. On the other hand, there seems to be clear indication that the scaling 
assumed in Eq. ([9]) holds in the bulk of the system. This scaling and the presence of 
gradients only in the direction perpendicular to the walls directly leads to the Fourier 
state discussed in Sec. [2j A different issue is to solve the Boltzmann equation to obtain 
the explicit expression for the one-particle velocity distribution of the state. The first 
Sonine approximation developed here clearly fails to give an accurate description of the 
tails of the velocity distribution. 

By construction, the Sonine approximation is expected to lead to much more 
accurate results for the low velocity moments than for the complete one-particle 
distribution function itself. Moreover these moments contain the most relevant physical 
information on the macroscopic properties of the system. A striking property of the 
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Figure 4. Steady dimensionless pressure, temperature, and density profiles in a system 
of hard disks with a = 0.99. The temperature has been scaled with some arbitrary 
reference value, Tr, actually the initial temperature of the system. The symbols are 
simulation data, the dashed straight line in the pressure profile is a guide for the 
eye, and the dashed line in the temperature profile is a linear fit of the data in the 
interval 50a < x < 700c The simulation parameters are: N = 3 x 10 3 , L x = 10 3 cr, 
Ly/ L x = 0.2. 



Fourier state investigated here is the relationship ( jT4l) between the temperature gradient, 
the pressure and the coefficient of restitution. Then, in Fig. [9j the MD simulation values 
obtained for the quantity I = 9 /pa in the bulk of the system has been plotted as a 
function of the restitution coefficient for the interval 0.9 < a < 1. These simulation 
values are seen to be in very good agreement with the theoretical prediction given by 
Eq. ( BSD , where the coefficients a i, &oi, and bio are given by the solution of the system of 
equations (1301 . (lA.lj) and (1A.2j) corresponding to the regular Fourier state. Also plotted 
in Fig. [9] is the quantity /, as computed from the hydrodynamic Navier-Stokes equations 
derived from the Boltzmann equation by means of the Chapman-Enskog procedure in 
the first Sonine approximation [2U EH] . These equations also admit a solution with the 
properties of the Fourier state [22J. The Navier-Stokes prediction is seen to be very close 
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Figure 5. The scaled distribution function for the same system as in Fig. [2l The 
symbols correspond to MD simulation data at three different values of the coordinate 
x, as indicated in the insert. The solid line is the theoretical prediction in the first 
Sonine approximation derived in the text, Eq. (|22|) . The dotted line is the Gaussian, 
included as a reference. 

to the result derived here. In fact, the difference between both predictions is smaller 
than the statistical uncertainties of the simulation data. 

Of course, another physically relevant moment is the heat flux q x , or its 
dimensionless form Q x defined in Eq. ( l24l) . This quantity is shown as a function of the 
restitution coefficient in Fig. [101 The observed qualitative behavior is easily understood, 
since it follows from Fig. [7] that the magnitude of the temperature gradient increases 
as a decreases. Again, there is a quite good agreement between the MD simulation 
results and the expression found in this paper, i.e. Eqs. fT24|) -f l26|) and 6 i an d b w 
obtained as indicated above and discussed in Sec. [3] for the regular solution. As in 
Fig. [9j the result obtained from the hydrodynamic Navier-Stokes equations, in the first 
Sonine approximation, for the Fourier state has also been included, and it differs very 
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Figure 6. The dimcnsionless heat flux Q x defined in Eq. (|24|) . as a function of the 
reduced coordinate x/a for the same system as in Fig. [TJ The symbols are MD 
simulation results and the dashed line indicates the theoretical prediction obtained 
in this paper. 



slightly from the solution derived here. 
6. Summary and discussion 

Two non-trivial solutions of the Boltzmann equation for smooth inelastic hard spheres 
or disks have been investigated. They describe a stationary state with constant 
pressure, temperature gradient in one direction, and no macroscopic mass flow. Quite 
peculiarly, the temperature profile is strictly linear and the heat flux can be expressed 
as proportional to it, without nonlinear contributions. In other words, the macroscopic 
state does not present any kind of rheological hydrodynamic effects. For this reason 
the state was termed Fourier state. Although no rigorous mathematical proof of the 
existence of the solutions is given, it has been shown that supposing they exist, a 
consistent expression for them can be found under well defined approximations. One of 
the solutions is singular in the sense that it does not reduce to any solution of the elastic 
Boltzmann equation when the limit of the coefficient of normal restitution going to unity 
is considered. It provides an example of a possible normal solution of the Boltzmann 
equation that is not captured by the Chapman-Enskog procedure to solve it. On the 
other hand, the other solution, termed regular, tends to an equilibrium Gaussian in the 
elastic limit. It is worth to stress that although both distribution functions are quite 
different, they lead to qualitatively similar macroscopic states. 

The singular state was never observed in the molecular dynamics simulations we 
carried out, while the theoretical predictions corresponding to the other solution are in 
good agreement with the simulation data, in the parameter region in which the Fourier 
state is accesible to the simulations, e.g. not too strong inelasticity. 

In this study, attention has been restricted to the first Sonine approximation, mainly 
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Figure 7. Steady dimensionless temperature and heat flux profiles in a system of hard 
disks with a = 0.9. The temperature Tr is the initial temperature of the system. The 
symbols are simulation data, the dashed line in the temperature profile is a linear fit 
of the data in the interval 250(7 — 550<r, and the dashed line in the heat flux profile is 
the theoretical prediction in this paper. The simulation parameters are: N — 3 x 10 3 , 
L x = 10 3 a, L y /L x = 2 . 



for the sake of simplicity and also in order to derive results in an analytical form. 
From this starting point, two extensions are possible. The first one is to incorporate 
more terms in the Sonine expansion of the distribution function, solving the resulting 
equations numerically. A way of implementing this is by following the method developed 
in [29] , that is tailored to solve the Boltzmann equation without considering any gradient 
expansion. The second extension refers to the tails of the distribution function. An 
asymptotic analysis of the inelastic Boltzmann equation for the Fourier state shows 
that for large values of the component of the velocity in the opposite direction to the 
increase of the temperature, the marginal velocity distribution exhibits an algebraic tail 
[T8] . This feature can also be analyzed in an numerically exact way by using the method 
presented in ref. [29]. This would allow, for instance, to predict the amplitude of the 
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Figure 8. The scaled distribution function for the same system as in Fig. The 
symbols correspond to MD simulation data at three different values of the coordinate 
x, as indicated in the insert. The solid line is the theoretical prediction in the first 
Soninc approximation derived in the text, Eq. (I22p . The dotted line is the Gaussian, 
included as a reference. 
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Figure 9. Dimensionless ratio I = 6/pa between the temperature gradient 9 and 
the pressure p in the Fourier state, as a function of the restitution coefficient a, for a 
system of inelastic hard disks. The symbols are from MD simulations, the solid line 
is the theoretical prediction derived in this paper, and the dashed line is the result 
obtained from the inelastic Navier-Stokes equations. 

algebraic decay. 

Future work should also attempt to put the existence of the Fourier state solution(s) 
of the Boltzmann equation on a more rigorous mathematical basis. Even if only the 
regular solution exists and this happens in some well defined limit, it should be an 
important result, since it would provide an explicit non-trivial normal solution of the 
Boltzmann equation. Moreover, the Fourier state also offers the opportunity of studying 
hydrodynamic fluctuations and correlations in a far from equilibrium state, without 
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Figure 10. Dimensionless heat flux Q x in the Fourier state as a function of the 
restitution coefficient, for a system of inelastic hard disks. The symbols are from MD 
simulations, the solid line is the theoretical prediction derived in this paper, and the 
dashed line is the result obtained from the inelastic Navier-Stokes equations. 

resorting to expansions in the gradients of the fields. 
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Appendix A. Two moment equations for the Fourier state in the first 
Sonine approximation 

Multiplication of Eq. ([TBI) by c z x and integration over c leads to 

64d(d + 2)(d + 4)(1 - a) [1 + 2(d - l)a 01 ] 

- (d - 1) [13 + 51a - 15(2 - 33ad - 4d 2 (l + 3a)] b 2 01 

- 9 [75 + 4d(21 + Ad - 9a) - 139a] b 2 l0 

- 6(d - 1) [29 + 8d(5 + d - 3a) - 93a] & i^io = 0. (A.l) 

In a similar way, multiplication by c 2 c x yields 

64d(d + 2)(d + 4)(1 - a) [d + 2 + (d + l)(d + 4)a 01 ] 
+ (d - 1) [262 - 390a - 7d + 231ad - d 2 (193 + 44d - 3(43 + 4d)a)] b 2 m 

- 9 [246 + 235d + 44d 2 - (22 + 3d) (17 + 4d)a] bj 

- 6(d - 1) [250 - 378a + 237d + 44d 2 - 3d(47 + 4d)a] b 01 b 10 = 0. (A.2) 

In the calculations leading to the above relations, Eq. (f29|) has been used and only 
terms up to second degree in aoi , &oi , and b w have been kept as discussed in Sec El The 
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circles plotted in Fig. [TJhave been obtained by solving Eqs. (1501) . (lA.ip . and flA.2[) for 
al = 2. 
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